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ABSTR  ACT(U) 

A  procedure  for  generating  disjoint  streams  of  pseudo-random  numbers  for  use 
in  discrete  event  simulation  is  presented.  The  procedure  facilitates  the 
synchronisation  of  variates  from  run  to  run  and  can  be  employed  in  conjunction 
with  several  variance  reduction  techniques.  It  is  based  on  Schrage's 
implementation  of  the  Lehmer  generator  and  includes  enhancements  to 
remove  the  restriction  on  usable  multiplier  values.  Seeds  are  efficiently 
computed  by  a  routine  that  exploits  the  number-theoretic  properties  of 
congruences.  Machine-portable  Fortran  routines  are  listed  for  use  on  machines 
capable  of  performing  31-bit  signed  integer  arithmetic. 
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1 .  INTRODUCTION 

Disjoint  streams  (viz  sequences)  of  pseudo- random  numbers  are  frequently 
employed  when  the  performances  of  several  systems,  system  configurations,  or 
operational  policies  are  compared  using  discrete  event  simulation.  Disjoint 
streams  are  employed  in  conjunction  with  several  variance  techniques, 
including  Common  Random  Numbers  and  Antithetic  Variates (ref . 1) ,  where  the 
objective  is  to  extract  in  the  minimum  possible  time  (or  run  length) 
confidence  interval  estimates  for  model  response  variables.  These  techniques 
require  a  generator  that  can  produce  independent  streams  and  reproduce  (viz 
synchronise)  them  from  run  to  run. 

During  the  course  of  building  a  GPSS/H  simulation  model  of  a  local  area 
necwork(ref . 10)  the  author  found  it  necessary  to  invoke  an  external  random 
number  generator  for  sampling  continuous  multi-parametric  distributions,  a 
facility  lacking  in  GPSS/H.  There  is  also  a  problem  with  the  high-order  serial 
correlation  properties  of  the  Tausworthe  random  number  generator  incorporated 
in  release  0.99  of  GPSS/H  (reference  1,  p  202). 

The  LCG,  or  Lehmer  generator (ref . 8) ,  is  the  traditional  source  of  pseudo¬ 
random  numbers.  The  prime  modulus  multiplicative  LCG  (PMMLCG) ,  described  by 
Hutchinson(ref  .4) ,  is  one  of  the  most  popular  forms  of  the  LCG  and  is  noted 
for  its  simplicity.  The  statistical  properties  of  the  PMMLCG  with  modulus 
231-1  are  the  also  best  known. 

The  PMMLCG  produces  a  sequence  of  non-negative  integers  <X1,X2, . >  from  a 

specified  initial  non-negative  integer,  X0 ,  using  the  recursive  linear 
congruence  formula 


X  , ,  =  (MX  )  mod  N,  for  n  =  0,1,... 

n+1  n  ’ 


(1) 


where  the  multiplier  (M)  and  the  modulus  (N)  are  integers  and  N  is  prime.  The 
symbol  "mod  N"  denotes  the  remainder  obtained  from  division  by  N  (see 
Appendix  I,  Definition  1).  Provided  M^N  and  X0*0  (mod  N) ,  the  values  of  the 

generated  integers  will  lie  in  the  closed  interval  [1,N-1].  The  Lehmer 
sequence  is  cyclical  with  a  period  dependent  on  M.  In  the  particular  case 
where  the  multiplier  is  a  primitive  root  of  the  modulus,  that  is,  if  the 

multiplier  satisfies  MN  1  =  1  (mod  N)  and  M^l  (mod  N)  for  0<t<N-l,  the  PMMLCG 
has  full  period  (N-l)  and  produces  (N-l)  unique  variates  in  the  cycle. 
A  stream  of  real  numbers  <Zt , Z2 ,....> ,  whose  values  are  asymptotically 

uniformly  distributed  in  the  open  interval  (0,1),  is  then  simply  obtained  by 
scaling,  viz  =  X^/N.  Variates  with  other  distributions  may  be  obtained  by 

transformation  (see  reference  1  pp  306  to  347  and  reference  7  pp  240  to  278). 

The  Lehmer  generator  is  featured  in  many  of  the  standard  mathematical 
subroutine  libraries  (eg  GGUBT  in  IMSL(ref . 15) ,  G05CAF  in  NAG(ref.l3)  and 
MTH$RANDOM  in  VAX/VMS(ref . 16) ) .  Recent  studies  by  Fishman  and  Moore(ref.2) 
suggest  that  the  statistical  properties  of  existing  PMMLCG  implementations  may 
not  be  optimal.  They  recommend  the  use  of  alternative  multipliers. 

An  efficient,  machine-portable  computer  program  capable  of  generating 
synchronised  disjoint  streams  of  pseudo-random  numbers  on  32-bit  machine 
architectures  is  described  in  this  paper.  It  is  based  on  the  well-known 
Schrage  algorithm(ref . 11) ,  a  Fortran  program  which  implements  the  PMMLCG  very 
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efficiently  on  32-bit  machines  by  using  integer  arithmetic.  Improvements  have 
been  made  to  increase  the  numerical  precision  of  the  calculations  to 
accommodate  the  multiplier  values  recommended  by  Fishman  and  Moore. 

The  remainder  of  this  paper  is  organised  as  follows.  In  Section  2  alternative 
ways  of  allocating  random  numbers  to  stochastic  processes  in  discrete  event 
simulation  models  are  illustrated  using  a  tandem  queuing  model.  Generation 
methods  based  on  interleaved  and  sequential  sampling  of  the  Lehmer  sequence 
are  discussed.  The  latter  method  forms  the  basis  of  the  two  algorithms 
subsequently  described  in  detail. 

The  first  algorithm,  described  in  Section  3,  is  the  improved  version  of 
Schrage's  random  number  generator.  The  second  algorithm,  described  in 
Section  4,  is  used  to  initialise  the  seeds  for  the  generator.  Provided  an 
upper  bound  can  be  specified  on  the  number  of  required  variates,  the 
initialisation  procedure  will  compute  seeds  that  enable  the  random  number 
generator  to  produce  disjoint  streams  in  accordance  with  the  sequential 
sampling  method. 

The  performance  of  each  algorithm  is  assessed  in  Section  5.  The  speed  of  the 
improved  Schrage  algorithm  is  compared  to  the  speed  of  the  original  Schrage 
algorithm  and  the  IMSL  GGUBT  algorithm,  and  the  time  taken  by  the 
initialisation  algorithm  to  compute  the  seeds  is  compared  to  the  projected 
time  to  compute  the  seeds  by  the  repeated  application  of  the  imp: jved  Schrage 
algorithm. 

This  research  was  undertaken  with  the  sponsorship  of  the  Director,  Naval 
Combat  Systems  Engineering,  Royal  Australian  Navy  (Task  No.  NAV87/226). 


2.  GENERATION  METHODS 

Variates  obtained  from  the  Lehmer  generator  can  be  assigned  to  stochastic 
processes  in  a  model  from  a  single  stream,  but  the  variates  may  be  correlated 
in  the  sense  that  if  a  process  is  modified  in  a  subsequent  simulation  run,  the 
sequence  of  variates  assigned  to  other  processes  may  be  altered.  Consequently 
this  method  cannot  be  used,  except  in  the  case  of  extremely  simple  models  (eg 
single  queues),  in  conjunction  with  the  most  popular  of  the  variance  reduction 
techniques  (viz  Common  Random  Numbers  and  Antithetic  Variates).  This  is 
illustrated  in  the  following  example. 

Consider  a  dual  tandem  queue.  Denote  the  arrival  of  customer  n  by  A^,  the 

commencement  of  service  at  service  centre  1  by  S  and  the  commencement  of 

n 

service  at  service  centre  2  by  Dn,  and  suppose  the  following  sequence  of 

events  occurs:  <Aj ,SX ,k2  ,Di ,S2 ,A3 ,D2 ,S3 ,A4 ,A5 ,S4 , . . .>.  If  numbers  in  the 

Lehmer  sequence  are  chronologically  assigned  to  processes  and  Zx  is  assigned 

to  the  first  event  (viz  A3),  then  customer  inter-arrival  times  and  centre  2 

service  times  would  be  computed  from  <Z3  ,Z3 ,Z6 ,Z9 ,Z10 , . . .>  and  <Z4,Z7,...> 

respectively.  Now,  consider  a  subsequent  replication  in  which  the  service 
times  of  centre  1  are  changed  and  the  resulting  event  sequence  is 
<A3  ,S  i  ,D i , A2  , S2  ,  A3  ,Dj  , S3 , A,,  ,S(, ,  As ,  .  .  .  > .  Customer  inter-arrival  and  centre  2 

service  times  would  be  computed  from  <Z 3 , Z4 , Z6 , Z9 , Zx  x  , . . . >  and  <Z3,Z7,...> 

respectively,  which  are  clearly  different  from  before. 

In  order  to  assure  that  processes  are  independent  from  run  to  run,  it  is 
necessary  to  allocate  a  disjoint  stream  to  each  process.  This  may  be  achieved 
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by  interleaved  sampling  or  by  sequential  sampling  of  the  Lehmer  sequence.  In 
the  previous  example  interleaved  sampling  involves  allocating  stream  1, 
comprising  <^i ,Z1+k>Z1+2^, . . .>,  to  the  arrival  process;  stream  2,  comprising 

<Z2 »^2+k’^2+2k’ ' ' ,>*  t0  t*ie  centre  1  service  process;  and  stream  3,  comprising 

<ZS  >Z^+jc>Z2+2jc»  •  •  •>»  to  the  centre  2  service  process.  Provided  k  is  greater 

than  or  equal  to  the  number  of  streams  (viz  k>3),  identical  random  arrival  and 
centre  2  service  times  may  be  assigned  to  the  respective  customers  in  each  run 
regardless  of  the  distribution  of  centre  1  service  times. 

The  main  problem  with  interleaved  sampling  is  that  a  variate  Z^+n^  in  stream  i 

cannot  be  very  efficiently  generated  from  its  progenitor  Z^+^n  ^  because 

between  one  and  a  (k-1)  intermediate  variates,  depending  upon  the  method  usei. 
has  to  be  generated  (and  wasted)  to  obtain  each  variate. 


Sequential  sampling  is  the  best  solution  for  generating  disjoint  streams 
because  only  one  application  of  the  PMMLCG  congruence  formula  (equation  (1)) 
is  needed  to  generate  each  variate.  In  this  case  variates  are  allocated  in 
contiguous  subsequences  as  follows.  Stream  1  comprises  <Zj ,Z2 ,Z2 , . . . ,Z .>, 

stream  2  comprises  <Z  ,Z  ,Z.  >  and  stream  3 

j+1  j+2  k 


comprises 


<Zk+l’Zk+2’ 


,Z  >.  Provided  not  more  than  j,  k-j  and  m-k  variates  are 

m  j  >  j 


respectively  required  from  streams  1,  2,  and  3  during  the  simulation  runs,  the 
streams  will  be  disjoint. 


The  initial  integer  variates,  or  seeds  (Z0, 


Z  .  and  Z.  )  , 
J  k 


must  be  computed  to 


ensure  that  the  streams  are  disjoint.  The  computations  need  only  be  performed 
once  prior  to  the  first  run.  The  computational  problem  is  similar  to  that 
encountered  previously  in  the  interleaved  sampling  method,  but  it  is  more 
constrained  because  it  is  infeasible  to  use  the  PMMLCG  formula  when  many 
thousands,  or  possibly  millions,  of  intermediate  variates  must  be  computed. 


3.  THE  IMPROVED  SCHRAGE  ALGORITHM 

The  naive  implementation  of  equation  (1)  using  P-bit  precision  positive 

p 

integers  will  incur  overflow  unless  the  product  M(N-l)  <2-1.  Many 
generators  of  interest  do  not  satisfy  this  criterion.  A  simple  solution  is  to 
employ  multiple-precision  floating  point  variables,  rather  than  integers, 
provided  the  precision  of  the  mantissa  is  adequate.  (On  a  given  machine, 
numbers  can  usually  be  stored  with  higher  precision  in  floating  point 
representation.)  However,  a  more  efficient  solution  may  be  achieved  using 
multiple-precision  integer  arithmetic  where  each  variate  and  the  constant 
multiplier  is  expressed  as  an  expansion  in  the  powers  of  an  integer  constant, 

r  =  2  ,  such  that  k<P.  The  modulo  N  product  may  then  be  obtained  by 
manipulating  the  expansion  coefficients. 

The  algorithm  described  by  Schrage  employs  radix  216  integer  arithmetic  to 
solve  equation  (1)  on  a  32-bit  machine  (viz  P  =  31)  using  a  31-bit  modulus, 
N  =  2J1-1,  and  a  15-bit  multiplier,  M  =  7s  =  16807.  Note  that  although  the 
Schrage  algorithm  employs  a  primitive  root  multiplier,  any  other  multiplier 
will  also  work,  provided  it  has  less  than  16-bit  precision  (viz  a  value  less 
than  32  768). 

Following  recent  studies  by  Fishman  and  Moore(ref.2)  there  has  been 
considerable  interest  in  the  use  of  larger  multipliers.  These  authors 
exhaustively  tested  the  statistical  properties  of  the  PMMLCG  with  N  =  2,l-l 
for  every  possible  primitive  root  multiplier  and  published  a  list  of  the  414 
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best  multipliers.  All  of  the  multipliers  exceed  23  bits  in  precision  and  thus 
cannot  be  employed  in  the  original  Schrage  algorithm.  Some  of  the  listed 
multipliers  in  current  usage  include  M  =  630,360,016  in  the  SIMSCRIPT  II 
simulation  language(ref .5) ,  M  =  397,204,094  in  the  IMSL  Fortran  subroutine 
library(ref . 15)  and  SAS  language(ref . 14) ,  and  M  =  742,938,285  in  the  GPSS/H 
simulation  language(ref . 3) . 

The  improved  Schrage  algorithm,  described  below,  is  similar  to  the  original 
algorithm  except  for  modifications  that  enable  the  use  of  31-bit  precision 
multipliers . 

Expanding  upon  the  original  Schrage  procedure,  the  integer  variable  (X^)  and 

the  multiplier  (M)  in  equation  (1)  are  represented  as  quadratic  and  linear 
polynomials  in  the  radix  r  =  215,  viz 


n 


=  U(2V  +  U(»t  4  U<0) 

n  n  n 


and 


M  =  v(1)r+v(0). 

n  n 

The  integers  u^\  an<*  Un^  ’  (corresponding  respectively  to  the 

15  low  order  bits,  the  15  middle  order  bits,  and  the  high  order  bit  of  X  )  and 

v<°>  and  v(1)  (corresponding  respectively  to  the  15  low  order  bits  and  the 


16  high  order  bits  of  M)  are  given  by: 


u 


(0)  _ 


u 


n 

(1) 

n 


(X  mod  r2 )  mod  r , 
n 


=  [(Xn  mod  r2) .  r  1 J, 


u 


(2)  _ 


n 

.(0) 


-  IV  *J. 


=  M  mod  r, 


and 


v 


(1) 


where  L-J  represents  the  largest  integer  smaller  than,  or  equal  to  (•)  (see 
Appendix  I,  Definition  2). 

The  product  (MX^)  may  also  be  expanded  as  a  polynomial  of  degree  four  in  r: 


MX 


„<4V  4„<3V  4„<2V  4„<1>r+u«» 

n  n  n  n  n 


(2) 
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where  the  coefficients  w^^ , . ,w^4\  corresPond  respectively  to  the 

15  low  order  bits,  the  three  15-bit  groupings  of  the  middle  order  bits,  and 
the  two  high  order  bits  of  the  62-bit  product.  Each  of  the  coefficients  can 
be  efficiently  computed  from  the  coefficients  of  the  integer 

variable  (X  )  and  the  multiplier  (M)  using  the  classical  multiplication 
algorithm  described  by  Knuth  (reference  6,  p  253). 

Simulated  division,  described  by  Payne  et  al(ref.9),  is  then  employed  to 
compute  the  product,  modulo  N.  Letting 

5n  =  L^n'  r"2/2J  =  Wn4)‘r2/2  +  Wn3)'r/2  +  LWn2)/2J  (3) 

from 

Y  if  Y  >0 
^  n  n 

Y  +  N  if  Y  <0 

n  n 


MX  -  N.(£  +1).  (4) 

n  n 


the  variable,  X  ,,,  is  obtained 
n+1 


n+1 


where 


Y  = 
n 


Substituting  equations  (2)  and  (3)  into  equation  (4),  and  collecting  terms  in 
the  powers  of  r,  Y^  is  then  also  obtained  in  terms  of  the  expansion 

coefficients  of  the  product,  viz 


Y 

n 


(w^4)  +  2w^2)  -  4(1  +  LWn2)/2J)) ' (r2/2) 
+  (w£3)  +  2w^1)).(r/2)  +  (1  +[_Wn2)/2J 


+  w 


(0) 


)• 


Note  that  each  of  the  bracketed  terras  in  the  above  expression  may  be  computed 
and  summed  without  overflow  because  their  values  lie  within  the  range  of 
31-bit  signed  integer  representation,  viz  [-23 l  +  l,231-l]  . 

The  Fortran  function  DRANX,  described  in  Appendix  II,  is  an  implementation  of 
the  improved  Schrage  algorithm.  Note  that  Knuth's  notation  has  been  retained 
for  the  expansion  coefficients  in  the  Fortran  program.  The  coefficients  are 
therefore  indexed  in  reverse  order,  commencing  at  one  instead  of  zero. 


4.  INITIALISATION  PROCEDURE 

The  number  of  iterations  of  equation  (1)  needed  to  obtain  a  seed  from  a 
reference  variate  X0,  in  the  Lehmer  cycle  is  called  the  offset,  0(s),  of 

stream  s  (s  =  1,2,..., S).  Assuming  the  offsets  are  totally  ordered  such  that 
O<0(1)<0(2)<. . .<0(S) ,  it  is  only  necessary  that  the  number  of  variates  in  any 
stream  s  does  not  exceed  (0(s+l)-0(s))  to  guarantee  that  the  streams  will  be 
disjoint  (viz  sequences  of  variates  do  not  overlap). 
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The  period  of  any  primitive  root  PMMLCG  with  modulus  N  =  231-1  should  provide 
a  sufficient  number  of  streams  of  the  required  length  for  practical  purposes. 
These  generators  can  be  configured  to  produce,  for  example,  more  than  a 
thousand  streams  with  one  million  variates  per  stream,  or  more  than  a  hundred 
streams  with  ten  million  variates  per  stream.  Generators  with  non-prime 
multipliers  can  also  be  employed  but,  because  their  cycle  lengths  are  shorter, 
larger  offsets  must  be  specified  to  obtain  the  same  number  of  variates  per 
stream. 

The  Lehmer  generator  itself  cannot  be  used  to  compute  seeds  because  it  is  far 
too  inefficient  to  compute  each  of  the  (8(s)-l)  intervening  variates. 

It  can  be  shown  by  induction  (Appendix  I,  Theorem  1)  that  equation  (1)  is 
equivalent  to 


Xn  =  (MnX0)  mod  N 


so  that  a  seed  is  given  by: 


X,(s)  =  (M0(S)XO)  mod  N 


(5) 


Equation  (5)  suggests  a  single  step  method  for  computing  seeds  but  it  is 
impractical  because  the  powers  of  large  multipliers  cannot  be  efficiently 
computed  or  stored. 

A  new  algorithm  which  exploits  the  number-theoretic  properties  of  congruences 
to  compute  the  seed  of  each  random  number  stream  is  described  below.  It 
involves  the  relatively  simple  computation  of  a  limited  number  of  intervening 
variates.  The  Fortran  subroutine  STREAM,  described  in  Appendix  II,  is  based  on 
this  algorithm. 

Firstly,  it  can  be  shown  (Appendix  I,  Theorem  2  and  Corollary  1)  that: 

X2n  =  (Q.(X^)  mod  N)  mod  N  ,  (6) 


where  the  factor  Q  is  the  integer  solution  of  the  Diophantine  equation: 


N.L  +  Q.X„  =  1 


(7) 


and  L  is  also  an  integer  solution.  Provided  that  N  and  X0  are  relatively 

prime  (ie,  the  greatest  common  divisor  of  N  and  X0  is  unity),  a  numerical 

solution  of  equation  (7)  can  be  obtained  using  an  extension  to  the  classical 
Euclid  algorithm  described  by  Knuth  (reference  6,  p  352).  Note  that  if  Q  is 
negative,  N  must  be  added  to  it  prior  to  substitution  into  equation  (6)  in 
order  to  maintain  X2n  non-negative. 

Next,  a  unique  sequence  A(k„)  can  be  defined  on  the  set  {0,1}  for  each 
integer,  k0 : 


a . 

=  <a.=l-(k.  mod  2):k.+1=l+(k.-2)/2  1 ,k.+1>l , i=0 , 1 1> 


A(k0 ) 


(8) 


7 


WSRL-TN-46/89 


The  elements  of  the  sequence  will  be  referred  to  as  (binary)  indicators .  The 
number  of  indicators,  £,  which  is  a  function  of  k0 ,  is  a  maximum  when  (k„+l) 

p 

is  a  second  or  higher  power  of  two;  ie,  if  k0  =  (2  -1),  then  £=2(P-1).  Thus, 

£  will  always  satisfy  £  <  2 ( log2 (k0+l ) - 1 )  provided  k0>l.  A  maximum  of 

60  indicators  is  required,  for  example,  to  represent  a  31-bit  integer  (viz 
P  =  31) . 

Now,  by  letting  ka  =  9(s),  equation  (5)  may  be  resolved  into  a  set  of  (£+1) 

equations.  Each  equation  is  a  congruence  relation  involving  an  indicator  of 
8(s),  as  follows: 

S£+* 

X  =  (M  XB)  mod  N  (9) 

k£ 


Xk  =  (M  £_1  Xk  )  mod  N  (10.1) 

£- 1  K£ 


a„+l 

X  (s)  =  X.  =  (M  X,  )  mod  N  (10.  £) 

U  K  p  K.  ^ 

where  a^  =  0  (and  k^  2  1),  and  k0  >  k2...>  k^.  Note  that  equation  (9)  may  be 
simplified  to 

X,  =  (MX0)  mod  N  =  X2,  (11) 

k£ 

where  Xx  is  the  first  variate  in  the  sequence  generated  by  equation  (1). 

Equations  (9)  to  (10. £-1)  represent  the  £  intervening  variates  that  must  be 
computed  between  the  reference  variate  and  the  seed. 

When  combined  with  equations  (1)  and  (6),  equations  (10.1)  to  (10.  £)  may  be 
expressed  by  a  single  recursion  formula, 

1-a.  a,  1+a . 

X.  =  (M  XQ  1(Xi+11)  mod  N)  mod  N  (12) 

where  i  =  £,£-l,...,l  and  =  Xt  and  XD(s)  =  Xj. 

The  seed  X0(s)  for  stream  s  offset  by  9(s)  is  therefore  obtained  from  the 
following  four  step  algorithm: 

(1)  compute  the  value  of  Q  in  equation  (7)  for  the  selected  PMMLCG  (and  if 
Q<0  replace  Q  by  Q+N) ; 

(2)  compute  =  Xx  from  equation  (11); 

(3)  determine  the  sequence  of  binary  indicators  <a0 ,ax  , . . . ,a^_ 
corresponding  to  the  integer  k0  =  9(s)  using  equation  (8),  and 
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(4)  iteratively  solve  equation  (12)  starting  with  i  =  l.  The  final  iterate 
(with  x  =  0)  is  the  required  seed. 


5.  COMPARATIVE  SPEED  OF  PMMLCG  ALGORITHMS 

The  mean  execution  times  of  the  Fortran  procedures  DRANX  (using  the  IMSL 
multiplier  value)  and  RAND,  which  implement  the  improved  and  original  Schrage 
random  number  generators  respectively,  were  measured  on  a  DEC  VAX-8200 
mainframe  computer  and  compared  to  the  execution  time  of  the  IMSL  Fortran 
random  number  generator  GGUBT(ref . 15) .  GGUBT  employs  double-precision  (64-bit 
mantissa)  floating  point  arithmetic  to  solve  equation  (1)  using 

X  ^  =  (fv(1).((rX  )  mod  N) )  mod  N  +  [v(2)X  ]  mod  N)  mod  N  . 
n+1  1  n  n 

The  advantage  of  using  integer  rather  than  floating  point  arithmetic  is 
evident  on  comparing  DRANX  and  GGUBT  execution  times.  The  penalty  for 
increasing  the  precision  of  the  multiplier  in  DRANX  is  reflected  in  the  RAND 
and  DRANX  execution  times.  For  each  call  DRANX  took  about  half  the  time 
(202  ±  3  ys ,  99%  confidence  interval  (c.i.))  of  GGUBT  (411  ±  28  ys ,  99%  c.i.), 
whereas  DRANX  took  about  twice  the  time  of  RAND  (119  ±  3  ys ,  99%  c.i.). 

The  total  execution  time  of  the  initialisation  procedure  STREAM  depends  upon 
the  number  of  streams  required  (S)  and  is  a  logarithmic  function  of  the 
individual  stream  offsets.  It  is  given  by  860(log2[  II  (8(s)+l)]-S)  (in  ys) 

l<s<S 

provided  8(s)£l  for  all  streams  s  =  1,2 . S.  The  constant  of 

proportionality  was  determined  by  regression  analysis  of  the  execution  times 
measured  on  the  VAX.  It  takes  into  account  that  steps  (1)  and  (2)  of  the 
algorithm  need  only  be  performed  once,  irrespective  of  the  number  of  streams. 
The  maximum  execution  time  occurs  with  an  offset  6(s)  =  N  when  the  number  of 
iterations  of  equation  (12)  is  a  maximum,  viz  i  =  60.  The  effectiveness  of 
the  algorithm  is  apparent  when  the  maximum  execution  time  (26  ms)  of  STREAM 
(with  S  =  1)  is  compared  with  the  estimated  time  taken  (133  h)  to  generate  the 
same  seed  from  equation  (1)  using  DRANX. 


6.  PORTABILITY 

The  programs  DRANX  and  STREAM  have  been  successfully  tested  on  a  variety  of 
computers,  operating  systems  and  Fortran  compilers.  Further  details  are 
provided  in  Appendix  II.  A  user  may  confirm  the  correct  operation  of  the 
procedures  by  running  the  test  program  provided  and  comparing  results. 
Although  the  programs  were  designed  for  use  on  32-bit  machines,  they  should 
work  on  any  16-bit  machine  (and  possibly  some  8-bit  machines)  provided  the 
Fortran  compiler  that  is  used  has  a  switch  for  IN'TEGER*4  format.  The  programs 
were  shown  to  work  on  the  16-bit  IBM  PC/AT.  Program  execution  times  were  not 
measured  on  the  PC/AT  but  they  are  expected  to  be  somewhat  longer  than  on  a 
32-bit  machine  operating  at  the  same  clock  speed  due  to  the  additional  machine 
instructions  needed  to  manipulate  and  store  integer  data.  The  execution  time 
of  DRANX  on  the  16-bit  machine  should  be  comparable  to  the  execution  time  of 
the  floating  point  GGUBT  routine  on  the  32-bit  machine. 


7 .  SUMMARY 

Disjoint  streams  of  pseudo-random  numbers  are  frequently  employed  when 
comparing  the  performances  of  several  systems,  system  configurations,  or 
operational  policies  using  discrete  event  simulation.  They  are  often  used  in 
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conjunction  with  variance  reduction  techniques  (eg  Common  Random  Variates  and 
Antithetic  Variates)  which  require  streams  to  be  reproduced  from  run  to  run 
(viz  synchronised  between  runs). 

The  allocation  of  disjoint  streams  to  stochastic  processes  in  a  simulation 
model  was  illustrated  for  a  simple  tandem  queue.  Two  methods  of  sampling  the 
linear  congruential  (or  Lehmer  generator)  sequence  to  obtain  reproducible 
streams  were  described. 

A  Lehmer  generator  algorithm  was  described  for  generating  run-synchronised 
disjoint  streams  on  a  computer  capable  of  performing  31-bit  signed  integer 
arithmetic.  The  Schrage  algorithm  was  selected  because  it  is  machine- 
independent  and  employs  very  efficient  integer  arithmetic.  Although  several 
prime-multiplicative  LCG  (PMMLCG)  algorithms  exist  in  mathematical  subroutine 
libraries,  none  of  them,  to  the  author's  knowledge,  employs  the  optimal 
multiplier  values  recommended  by  Fishman  and  Moore.  Modifications  to  the 
Schrage  algorithm,  which  increase  the  multiplier  precision  from  24  to  31  bits, 
facilitate  the  use  of  the  recommended  values. 

The  algorithm  was  described  mathematically  and  presented  as  a  Fortran  function 
(DRANX) .  The  real  variates  derived  from  this  function  are  asymptotically 
uniformly  distributed  in  the  unit  interval  (0,1)  and  can  be  transformed  to  any 
other  common  distribution. 

The  difficulties  involved  in  computing  the  appropriate  seeds  for  this 
generator  were  discussed.  An  efficient  iterative  procedure  that  exploits  the 
number-theoretic  properties  of  congruences  to  compute  the  seeds  was  described. 
It  entails  the  conversion  of  a  vector  of  stream  offsets  (denoting  the 
positions  of  the  seeds  relative  to  a  reference  point  in  the  Lehmer  cycle)  to 
seed  values.  This  procedure  was  implemented  as  a  Fortran  subroutine  (STREAM) 
and  is  used  to  initialise  the  random  number  generator.  Provided  the  user- 
specified  stream  offsets  satisfy  a  simple  ordering  relation,  the  streams  can 
be  guaranteed  to  be  disjoint  and  reproducible  between  runs. 

The  execution  times  of  the  improved  Schrage  algorithm  were  measured  and 
compared  with  the  original  Schrage  and  the  IMSL  GGUBT  algorithms.  Each  call  to 
GGUBT  took  about  twice  as  long  as  a  call  to  DRANX  since  GGUBT  employs  floating 
point  arithmetic,  whereas  DRANX  employs  integer  arithmetic.  DRANX,  on  the 
other  hand,  took  about  twice  as  long  as  the  original  Schrage  algorithm  due  to 
the  additional  complexity  in  using  31-bit  precision  multipliers.  The  maximum 
execution  time  of  STREAM  was  measured.  On  the  mainframe  test  computer  STREAM 
took  26  ms  against  an  estimated  133  h  CPU  time  needed  to  generate  a  seed  from 
a  reference  variate  by  repeated  calls  to  DRANX.  This  result  attests  to  the 
efficiency  of  the  initialisation  procedure.  The  portability  of  DRANX  and 
STREAM  have  been  verified  on  several  machines  (including  a  16-bit  machine), 
operating  systems  and  Fortran  compilers. 
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APPENDIX  I 

NUMBER-THEORETIC  DEFINITIONS  AND  PROOFS 


Definition  1 

Two  integers  a  and  b  are  congruent  modulo  m  if  their  difference  is  a  multiple 
of  a  non-zero  integer  m  (m  is  called  the  modulus)  viz.,  (a-b)  =  nm  for  some 

integer,  n.  The  congruence  relation  is  expressed  by  the  notation 

a  =  b  (mod  m) . 

Definition  2 

L*J  denotes  the  largest  integer  smaller  than,  or  equal  to  a  real  parameter  (•) 
(viz,  if  z  =  1+6,  where  I  is  an  integer  and  0£6<1,  then  LZJ  =  l)- 

Lemma  1 

If  a  and  m  (m^O)  are  integers,  then  a  mod  m  =  a-mj^a/mj. 

Proof:  Let  b  =  a-m  [a/mj .  Clearly  the  difference  of  a  and  b  is  an  integer 
multiple  of  m,  hence  b  =  a  (mod  m)  =  a-m|a/mj.  QED. 

Lemma  2 

If  the  integers  a  and  b  are  smaller  than  m,  and  x  =  (ab)  mod  m  and 

y  =  ((a  mod  m) . (b  mod  m))  mod  m,  then  x  =  y. 


).(b-m[b/ 

=  (ab-am[b/mj -bm  Ja/mJ  +m2  [a/mj  .  [b/mj  )  mod  m 

Since  a,b  <  m,  then  J^a/m  J  =  [b/mj  =  0  (by  Definition  2). 

Hence  y  =  (ab)  mod  ra  =  x.  QED. 

Theorem  1 

The  n-th  term  in  the  integer  sequence  <^j:J  =  generated  by 

equation  (1)  Xn  =  (MnX0)  mod  N.  This  is  a  well-known  result  (ref .  7 ,  p  222) 
which  can  be  proved  by  induction. 

Proof:  Xx  =  (MX0)  mod  N 

X2  =  (MX!)  mod  N  =  (M(MX0)  mod  N)  mod  N 

2 

=  (M  X,)  mod  N  (by  Lemma  2) 


mj )  mod  m  (by  Lemma  1) 


Proof:  y  =  ((a-m|a/mj 


Assume  X  =  (MnX#)  mod  N 
n 
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Consider 

Xn+1  =  (MXn)  mod  N  =  (M(MnX0)  mod  N)  mod  N 
=  (Mn+^X0)  mod  N  (again  by  Lemma  2) 

Hence,  by  induction, 

Xn  =  (MnX0)  mod  N  for  all  n  =  1,2, _  QED. 

Lemma  3 

If  X^  and  X^n  are  the  n-th  and  2n-th  terms  of  the  integer  sequence 

2 

<Xj :  j=l ,  2  ,....> ,  generated  by  equation  (1),  then  (XoX2n)  mod  N  =  (X^)  mod  N. 
Proof:  ^°^2n^  mo<*  ^  =  (Xo(M^nXo)  mod  N)  mod  N  (by  Theorem  1) 

=  ((MnXfl)  mod  N.(MnX0)  mod  N)  mod  N  (by  Lemma  2) 

=  (X*  )  mod  N.  QED. 

Theorem  2 

If  X^  and  X^n  are  terms  of  the  integer  sequence  <X . : j=l ,2 , . . . > ,  generated  by 
equation  (1),  then  X^n  may  be  expressed  as  a  function  of  X0  and  X^  by 
2 

X2n  ~  (Q(X  )  raoc*  mod  ^ >  where  Q  is  the  solution  of  the  congruence  relation 
(QX0 )  mod  N  =  1. 

Proof:  Suppose  X^n  =  (Q(X^)  mod  N)  mod  N 

then  (XQX  )  mod  N  =  (X0(Q(X*)  mod  N)  mod  N 

=  ((QX0)  mod  N.(X^)  mod  N)  mod  N  (by  Lemma  2) 

=  (X^)  moc*  ^  (hy  Lemma  3) 

The  supposition  is  correct  iff  (QX0)  mod  N  =  1.  QED. 
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Corollary  Is  The  solution  Q  of  (QX„)  mod  N  =  1  is  given  by  the  solution  of 
the  linear  Diophantine  equation:  NL  +  QX0  =  1,  where  Q  and  L  are  integers. 

Proof:  By  Definition  1,  QX0  -  1  =  -NL,  for  any  integer  L. 


Hence  NL  +  QX0 


1.  QED. 
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APPENDIX  II 
FORTRAN  ROUTINES 

In  Sections  II. 1  and  II. 2  two  ANSI  Fortran  66/77  compatible  procedures  (STREAM 
and  DRANX) ,  which  implement  the  multiple  stream  pseudo-random  number 
generator,  are  described.  Both  algorithms  are  designed  to  be  machine-portable 
and  suitable  for  computers  that  can  perform  31-bit  signed  integer  arithmetic. 
The  portability  of  the  routines  has  been  demonstrated  using  four  different 
compilers  and  three  machines/operating  systems  as  follows: 


Machine 

Operating  system 

Compiler 

IBM  3081 

TSO/E  2.0 
+  MVS/XA  2.1.2 

IBM  Fortran/Gl  l.l1 

DEC  VAX  8200 

VMS  4.7 

VAX  Fortran  4.8-276*  3 

IBM  PC/AT 

DOS  3.1 

Ryan-McFar land  Fortran  1.01 

II 

II 

Ryan-McFarland  Fortran  2.43 

1  ANSI  66  2  ANSI  66  using  /NOF77  compiler  option  3  ANSI  77 


Disjoint  streams  are  obtained  by  extracting  the  variates  in  each  stream  from 
different  segments  in  the  cycle  of  a  multiplicative  prime  modulus  linear 
congruential  (Lehmer)  generator  with  modulus  231-1.  The  Fortran  function 
DRANX,  which  generates  the  variates,  is  an  improved  version  of  Schrage's 
algorithm(ref . 12) .  The  numerical  precision  of  the  Lehmer  multiplier  has  been 
increased  to  31  bits  to  facilitate  the  use  of  the  optimal  multipliers 
recommended  by  Fishman  and  Moore(ref . 2) .  Prior  to  calling  the  generator,  the 
seeds  for  each  stream  are  initialised  by  an  efficient  procedure  (STREAM)  that 
exploits  the  number-theoretic  properties  of  congruences.  The  variates 
obtained  from  the  generator  are  asymptotically  uniformly  distributed  in  the 
open  interval  (0,1). 

A  FORTRAN  procedure,  listed  in  Section  II. 3,  demonstrates  the  method  of 
invoking  the  generator.  The  test  results  are  illustrated  in  Section  II. 4. 

II. 1  INITIALISATION  PROCEDURE  (STREAM) 

C 

C  SUBROUTINE  STREAM  IS  CALLED  ONCE  PRIOR  TO  INVOKING  THE  FUNCTION 
C  DRANX.  STREAM  OFFSETS  ARE  INPUT  IN  VECTOR  IVEC  WHICH  ARE 
C  CONVERTED  TO  SEEDS  UPON  RETURN. 

C 

C  ARGUMENT  DESCRIPTIONS: 

C 

C  NS  -  THE  MAXIMUM  NUMBER  OF  STREAMS  TO  BE  ESTABLISHED  (INTEGER) 

C 

C  IVEC  -  VECTOR  (INTEGER)  OF  DIMENSION  NS  USED  TO  INPUT  STREAM 

C  OFFSETS  AND  SUBSEQUENTLY  TO  STORE  THE  LAST  VARIATES 

C  GENERATED  IN  EACH  STREAM. 

C 

C  IREF  -  THE  REFERENCE  VARIATE  (INTEGER) 

C 

C  MULT  -  LEHMER  MULTIPLIER  (INTEGER) 

C 

C  COMMON  BLOCKS:  LABELLED  COMMON  BLOCK  /MPARM/  USED  INTERNALLY 
C  TO  PASS  DATA  FROM  STREAM  TO  FUNCTION  DRANX. 

C 

C  USAGE : 

C 
C 


/MPARM/ V1,V2 


noon 
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C 

C  VARIABLES: 

C 

C  VI  -  THE  REMAINDER  (INTEGER)  FROM  THE  DIVISION 

C  OF  MULT  BY  2**15 

C 

C  V2  -  THE  QUOTIENT  (INTEGER)  FROM  THE  DIVISION 

C  OF  MULT  BY  2**15 

C 
C 

C  OTHER  ROUTINES  CALLED:  IQFAC  AND  LMOD 
C 

SUBROUTINE  STREAM (NS , IVEC.MULT, IREF) 

IMPLICIT  INTEGER  (A-Z) 

DIMENSION  IVEC(NS) ,A(60) 

COMMON  /MPARM/V1 , V2 
DATA  B15/32768/ 

C  . .  DETERMINE  EXPANSION  COEFFS  OF  MULT  FOR  FUNCTION  DRANX 
Vl=MOD(MULT,B15) 

V2=MULT/B15 

C  ..  COMPUTE  Q  (STEP  1) 

Q=IQFAC(IREF) 

C  . .  COMPUTE  XI  (STEP  2) 

Xl=LMOD ( IREF , MULT) 

DO  10  S=1,NS 

C  ..  DETERMINE  INDICATORS  A(1),..,A(L)  FOR  STREAM  S  (STEP  3) 
KO=IVEC(S) 

IF  (KO.EQ.O)  IVEC(S)=IREF 
1=1 

20  IF  (KO.LE.l)  GOTO  30 

A(I)=1-MOD(KO,2) 

K0=l+(K0-2)/2**A(I) 

1=1+1 
GOTO  20 
30  1=1-1 

C  ..  SOLVE  EQUATION  (18)  FOR  VARIATES  (STEP  4) 

XHAT=X1 

40  IF  (I.LT.l)  GOTO  50 

XHAT= LMOD (MULT** ( 1 - A ( I ) )*Q**A ( I ) , 

X  LMOD ( XHAT , XHAT**A (I))) 

1=1-1 
GOTO  40 

50  IF  (KO.NE.O)  IVEC(S)=XHAT 

10  CONTINUE 
RETURN 
END 

FUNCTION  IQFAC  COMPUTES  THE  FACTOR  Q  IN  EQ.(13)  USING 
EUCLID’S  ALGORITHM  (USING  KNUTH'S  NOTATION,  SEE  KNUTH  P.352) 

FUNCTION  IQFAC (I REF) 

IMPLICIT  INTEGER  (A-Z) 

Ul=l 
U2=0 

U3=2147483647 
V1=0 
V2=l 
V3=IREF 
10  IF  (V3.EQ.0)  GOTO  20 
Q=U3/V3 
T1=U1-V1*Q 
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T2=U2-V2*Q 
T3=U3-V3*Q 
U1=V1 
U2=V2 
U3=V3 
V1=T1 
V2=T2 
V3=T3 
GOTO  10 
20  1QFAC=U2 

IF  (U2 . LT. 0) IQFAC=IQFAC+U 
RETURN 
END 
C 

C  FUNCTION  LMOD  COMPUTES  (IX*IY)  MOD  (2**31-1),  WHERE 
C  IX, IY< (2**31- 1)  USING  CLASSICAL  MULTIPLICATION  ALGORITHM  (USING 

C  KNUTH'S  NOTATION,  SEE  KNUTH  P.253)  AND  SIMULATED  DIVISION. 

C 

FUNCTION  LMOD (IX, IY) 

IMPLICIT  INTEGER  (A-Z) 

DATA  B14/16384/ .B15/32768/ .B29/536870912/ , B30/ 1073741824/ 

DATA  P/2147483647/ 

IXS=IX 

U1=IXS/B30 

IXS=IXS-U1*B30 

U2=IXS/B15 

U3=IXS-U2*B15 

V1=IY/B15 

V2=IY-V1*B15 

T=U3*V2 

K=T/B15 

W5=T-K*B15 

T=U2*V2+K 

K=T/B15 

W4=T-K*B15 

T=U1*V2+K 

K=T/B15 

W3=T-K*B15 

W2=K 

T=U3*V1+W4 

K=T/B15 

W4=T-K*B15 

T=U2*V1+W3+K 

K=T/B15 

W3=T-K*B15 

T=U1*V1+W2+K 

K=T/B15 

W2=T-K*B15 

W1=K 

TX=l+W3/2 

LM0D=(W5+TX)+B 14* (W2+2*W4 )+B29* (W 1+2*W3-4*TX) 

IF  (LMOD.LT.O)  LM0D=LM0D+P 

RETURN 

END 


n  n 
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I I. 2  RANDOM  NUMBER  GENERATOR  (DRANX) 

C 

C  FUNCTION  DRANX  RETURNS  A  DOUBLE  PRECISION  VARIATE  UNIFORMLY 

C  DISTRIBUTED  IN  THE  OPEN  INTERVAL  (0,1) 

C 

C  ARGUMENT  DESCRIPTION: 

C 

C 

C  I VAR  -  INTEGER  VARIATE  (ALTERED  BY  THE  ROUTINE) 

C 

C  COMMON  BLOCKS:  LABELLED  COMMON  BLOCK  /MPARM/  USED  INTERNALLY 
C  TO  PASS  DATA  FROM  SUBROUTINE  STREAM  TO  DRANX  (FOR 

C  DETAILED  INFORMATION  REFER  TO  SUBROUTINE  STREAM) . 

C 

C  OTHER  ROUTINES  CALLED:  NONE 

C 

FUNCTION  DRANX ( I VAR) 

IMPLICIT  INTEGER  (A-C.E-X) 

COMMON  /MPARM/V1 , V2 

C  ..  INITIALISE  B 14=2**14 ,  B15=2**15,  B29=2**29 ,  B30=2**30,  P=2**31-l 
DATA  B14/ 16384/, B15/32768/.B29/5368709 12/ ,B30/ 1073741824/ 

DATA  P/2147483647/ 

PASS  PARAMETERS  FROM  INITIALISING  ROUTINE 
COMPUTE  EXPANSION  COEFFICIENTS  OF  IVAR 
U1=IVAR/B30 
IVAR=IVAR-U1*B30 
U2=IVAR/B15 
U3=IVAR-U2*B15 

C  ..  COMPUTE  PRODUCT  OF  MULT  AND  IVAR  (KNUTH'S  NOTATION) 

T=U3*V1 

K=T/B15 

W5=T-K*B15 

T=U2*V1+K 

K=T/B15 

W4=T-K*B15 

T=U1*V1+K 

K=T/B15 

W3=T-K*B15 

W2=K 

T=U3*V2+W4 

K=T/B15 

W4=T-K*B15 

T=U  2*V  2+W3+K 

K=T/B15 

W3=T-K*B15 

T=l'l*V2+W2+K 

K=T/B15 

W2=T-K*B15 

W1=K 

TX=l+W3/2 

IVAR=(W5+TX)+B14*(W2+2*W4)+B29*(W1+2*W3-4*TX) 

IF  (IVAR.LT.O)  IVAR=IVAR+P 
C  . .  NORMALISE  VARIATE 

DRANX=IVAR*4 . 6566128752457969D- 10 

RETURN 

END 


oook-*  ui  non  noo 
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I I. 3  TEST  PROGRAM 

THIS  PROGRAM  DEMONSTRATES  THE  USAGE  OF  THE  ROUTINES  DRANX 
AND  STREAM.  THE  OUTPUT  LISTS: 

(1)  THE  SEEDS  COMPUTED  BY  THE  SUBROUTINE  STREAM  FROM  THE 
OFFSETS  ASSIGNED  TO  STREAMS  ONE  TO  FOUR  AS  FOLLOWS: 

STREAM  OFFSET 

1  99,999 

2  999,999 

3  9,999,999 

4  99,999,999 

(2)  THE  FIRST  10  SCALED  (0,1)  VARIATES  IN  EACH  OF  THE  STREAMS 
THAT  ARE  COMPUTED  BY  FUNCTION  DRANX. 

CONSTANT  PARAMETERS: 

MULT  -  LEHMER  MULTIPLER  (GGUBS  VALUE  397,204,094). 

IREF  -  REFERENCE  VARIATE  (GGUBS  SEED  1,891,356,973). 
IVEC(l) , . . . ,IVEC(4)  -  STREAM  OFFSETS  (SEE  NOTE  1  ABOVE) 

INTEGER  MULT , IREF , IVEC (4) , NS, S, LOOP 
DATA  MULT/397204094/ , IREF/1891356973/ ,NS/4/ 

DATA  IVEC/99999, 999999, 9999999, 99999999/ 

EXTERNAL  DRANX, STREAM 

. .  INITIALISE  4  STREAMS 

CALL  STREAM (NS, IVEC, MULT, IREF) 

. .  PRINT  SEEDS 

WRITE  (6,5) 

FORMAT (’  STREAM  SEED') 

WRITE  (6 , 10) (S , IVEC(S) , S=1 ,NS) 

0  F0RMAT(4X, I 1 , 7X, I 10) 

. .  PRINT  VARIATES 

WRITE  (6,15) 

15  FORMAT (1 OX,/ '  RANDOM  NUMBERS  GENERATED  BY  DRANX’// 

X'  POSITION  STREAM  1  STREAM  2  STREAM  3  STREAM  4') 

DO  30  LOOP=1 , 10 

WR ITE ( 6 , 20 ) LOOP , (DRANX ( I VEC ( S ) ) , S= 1 , NS ) 

20  FORMAT (3X, 12 , 6X,4(F9 . 7 , 3X) ) 

30  CONTINUE 

END 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 


21 


WSRL-TN-46/89 


I I. 4  TEST  RESULTS 

STREAM  SEED 

1  1496156467 

2  1851398463 

3  81679943 

4  287046138 


RANDOM  NUMBERS  GENERATED  BY  DRANX 


POSITION 

STREAM  1 

STREAM  2 

STREAM  3 

STREAM  4 

1 

0.9922353 

0.4811183 

0.1602005 

0.3419817 

2 

0.4851381 

0.1301162 

0.9993262 

0.3450323 

3 

0.2193673 

0.7638627 

0.5313926 

0.9414775 

4 

0.8110999 

0.3543463 

0.5686302 

0.3577577 

5 

0.8778924 

0.1577095 

0.0287373 

0.7472422 

6 

0.1481317 

0.5516356 

0.6789374 

0.8971976 

7 

0.3538447 

0.6098807 

0.9130332 

0.6669795 

8 

0.5466737 

0.4177297 

0.6134635 

0.6418861 

9 

0.3290880 

0.5880992 

0.2224111 

0.9184812 

10 

0.6634960 

0.3210178 

0.1404637 

0.5570008 
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